Boundary normal pressure-based electrical conductivity reconstruction for magneto-acoustic tomography with magnetic induction
1. IntroductionMagneto-acoustic tomography with magnetic induction (MAT-MI), first proposed by Xu and He in 2005,[1] is a noninvasive multi-physics imaging approach, which integrates the good contrast of electrical impedance tomography (EIT)[2–4] with the improved spatial resolution (mm-level) of ultrasound imaging.[5,6] MAT-MI was demonstrated to have the capability of electrical impedance contrast imaging[7,8] of biological tissues with conductivity differences. In MAT-MI,[9–16] a time-varying magnetic stimulation is used to induce eddy current inside conductive object. In the presence of a static magnetic field, mechanical vibrations are generated by the Lorentz force acting on the induced eddy current to produce detectable ultrasound. Through acoustic transmission in media, the magneto-acoustic signals can be acquired by the transducers placed around, which can be used to reconstruct two-dimensional (2D) images relating to conductivity contrast of scanned layer. Therefore, by combining magnetic fields with acoustic field,MAT-MI is often referred to as a hybrid imaging technology, which shows good potential applications in early diagnosis of tumor.
In the past decades, several numerical and experimental investigations were conducted to demonstrate the feasibility and performance of MAT-MI. In 2008, based on the theories of eddy current induction and acoustic transmission, Ma and He[17] reported that the acoustic waveform was generated at the conductivity boundary and the polarity of wave cluster reflected the direction of the exerted Lorentz force caused by conductivity variation along the normal direction of the induced eddy current. Hu et al.[18] verified that MAT-MI could be used to image tumor inhuman liver tissue in vitro. Then based on the radiation theory of acoustic dipole source, a tomographic conductivity reconstruction algorithm of MAT-MI was presented by Sun et al.[19] and also verified by experimental measurements. Both the configuration and the inner conductivity of the scanned layer of the object were reconstructed with improved spatial resolution. Mariappan and He[20] combined the vector source reconstruction algorithm with the relationship obtained using Helmboltz decomposition of the induced eddy current to estimate the conductivity distribution.Guo et al.[21] reported that the relative conductivity distribution could be mapped with reciprocal theorem in magneto-acousto-electrical tomography with magnetic induction (MAET-MI) modality. In 2016, Yu[22] proved by using an in vivo mouse model that high-frequency MAT-MI was able to resolve the boundary and the internal structures of tumor.
In previous studies, although the capability of distinguishing tissues with conductivity difference was demonstrated by the 2D or three-dimensional (3D) image reconstructions using the waveforms collected around the object, several problems existing in physical meaning and signal processing still need to be solved. (i) The physical meanings of acoustic pressure and waveform are confused in image reconstruction and the influences of magnetic fields and the transducer impulse response are often neglected. (ii) The physical meaning of pressure derivative in traditional MAT-MI[23,24] is not explained clearly and the conductivity reconstruction process is just mathematically computed but not analyzed in physics. (iii) The dimension of the reconstructed image is not the electrical conductivity, even for the normalized results. Whereas, it was reported by Zhou et al.[25] that acoustic pressure with sharp and clear boundary peaks could be detected by the transducer with strong directivity and wide bandwidth, reflecting the gradient of the induced Lorentz force accurately. Corresponding to the acoustic sources at conductivity boundaries, the boundary pressures along the normal direction (boundary normal pressure, BNP) could be retrieved through the analyses of amplitude and phase of collected wave clusters, and then the conductivity distribution of the scanned layer could be reconstructed accurately.
In this study, based on the analyses of magnetic excitation and acoustic vibration in conductive object, acoustic radiation in diffraction mode, and acoustic detection with a planar piston transducer, the inverse problem solution for acoustic source reconstruction is derived in 3D configuration with explicit formulae. By considering the strong directivity of the transducer, the physical meanings of detected acoustic pressure and its derivative are analyzed, proving that BNP could be extracted exactly by the abrupt change of pressure derivative. Then, based on the analysis of wave cluster, a simplified BNP-based electrical conductivity reconstruction algorithm is developed by using Hilbert transform and the clear physical meaning of pressure derivative relating to the property of conductivity variation is demonstrated. The proposed algorithm is also verified by numerical studies for a 2-layer eccentric cylindrical model. The simulations consist well with the experimental results in our previous studies,[14,15,19] and the reconstructed images are well consistent with the cross-sectional conductivity distributions of the model. The favorable results prove the feasibility of the proposed conductivity reconstruction algorithm based on BNP and suggest a new simplified imaging method for potential applications in MAT-MI.
2. Principle and methodThe schematic diagram of MAT-MI system is shown in Fig. 1. A conductive object is placed in a complex magnetic field along the z direction, which is composed of a static magnetic field 
and a pulsed magnetic field 
. The magnetic field is assumed to homogenously cover the entire object. Excited by 
, the eddy current 
can be induced at 
inside the object as closed loops in the xy plane. Then, the acoustic pressure 
can be generated by the induced Lorentz force 
acting on the eddy current. Due to the interaction of the eddy current with the static magnetic field, the medium vibration and acoustic transmission are produced. After propagating through the object and surrounding water, the pressure 
can be detected by the transducer at 
to map the electrical conductivity distribution of the scanned layer.
In order to simplify the theoretical derivation, the pulsed magnetic field is set to be a Dirac delta function 
, the induced eddy current can be described by 
for 
. By considering the induced Lorentz force 
, 
satisfies the wave equation[25]
where
and
are the Laplacian operator and the divergence operator,
is the convolution operator, and
c is the acoustic velocity. As reported in previous studies, the Lorentz force induced one-dimensional (1D) acoustic vibration can be treated as a dipole source and the radiation pressure can be described by
, where
θ is the radiation angel between the normal direction
of the eddy current and the transmission path
. Meanwhile, based on the reciprocal theorem for the transducer, the reception pattern
[26] should be considered, where
is the regular cylindrical Bessel function of first order,
a is the radius of the planar piston transducer,
is the wave number for the angular frequency
ω,
ka is the key parameter describing the directionality of the transducer, and
β is the angle between the normal direction
of the transducer surface and the transmission path
. By employing Green’s function, the acoustic pressure
detected at
can be calculated from
[27]
where
ξ is the pressure conversion efficiency of the transducer,
is the time delay for acoustic transmission from
to
. Finally, by considering the waveforms of the pulsed magnetic field
and impulse response
of the transducer, the collected acoustic waveform can be calculated by
.
[27]By defining
as the equivalent acoustic source, the inverse problem solution for MAT-MI can be solved in 3D diffraction mode. In Fourier transform, we obtain
,
and
for
and
. Then the Fourier transforms of
and
on variable
can be obtained. Further taking Fourier transform
,
[28] the detected acoustic pressure in the frequency domain can be written as
In cylindrical coordinates, 
can be expanded to[28]
and equation (
3) can be transformed into
where
,
and
represent the positions
and
in the cylindrical coordinate, respectively.
is the
m-th order Bessel function of the first kind and
is the
m-th order Hankel function.
[29,30]In order to extract 
from the right-hand side in Eq. (4), mathematical transformations are conducted according to the following sequence. First multiplying both sides of Eq. (4) by 
and integrating them over 
from 0 to 
, and then replacing 
with 
, we obtain
where
can be simplified into
only for
. Further multiplying both sides of Eq. (
5) by
and summing
n from
to
,we obtain
and it can be simplified into
Then, by multiplying both sides of Eq. (7) with 
and 
, and integrating them over μ from 0 to 
and 
from 
to 
, the resulting equation can be rewritten as
After integral transformation for the right-hand side, equation (8) can be simplified into
So the equivalent acoustic source 
inside the object can be restored from 
. For a cylindrical model without considering the z direction, we obtain 
and 
. By replacing 
with 
, and μ with k for 2D cylindrical measurement configuration, equation (9) is changed into
and then, the equivalent source inside the object is achieved as follows:
 | |
For biological tissues with low conductivity, 
[31] and only the acoustic sources generated at conductivity boundaries should be considered. 
can be calculated to be 
by the simplified gradient of the induced eddy current along the normal direction. For a practical transducer with strong directivity, only acoustic signals transmitted along the normal direction can be detected in a narrow angle scope, and 
can be simplified into 1 for 
. Meanwhile, for two points 
and 
as shown in Fig. 1, when the direction of the induced Lorentz force is along the normal direction of the transducer, the radiation and reception angles θ and β are both zero, denoting that the pressure detected at 
transmitted from 
is the BNP of the source with a maximal amplitude. Therefore, in 
, vertical pressure jumps are generated at the conductivity boundaries in the direction of the induced Lorentz force along the normal direction of the transducer for 
, representing the corresponding BNPs. With the derivative of pressure 
with respect to time t, the amplitudes and directions of vertical pressure jumps can be retrieved to restore BNPs generated at the corresponding conductivity boundaries. So equation (11) can be transformed into
 | |
When 
, the normal direction 
of the induced eddy current points to the normal direction of the transducer. Since the direction of BNP is along the radiation direction (
of the dipole source, 
can be changed into 
for the BNP transmitted from 
to 
along 
, where 
. So the integral over 
can be changed into ct, and the conductivity distribution can be reconstructed by[16]
where
Rm and
are the scanning radius and measured angle of the transducer. For a medium with the area and perimeter of
S and
L, the magnetic induction satisfies
and the reconstructed conductivity can be corrected into
[16]
In practical applications for MAT-MI, because the experimental waveform 
(system transfer function) is not accurate enough, the exact distribution of 
is difficult to de-convolute from 
, even for the widely used Wiener inverse filter[32] in ultrasound imaging. Whereas, by applying Hilbert transform,[33] BNPs at conductivity boundaries can be retrieved with accurate positions and amplitudes as well as vibration polarities from the wave clusters in 
, and then 
can be restored for further conductivity image reconstruction based on Eq. (14).
3. Numerical studiesIn order to verify the performance of the proposed BNP based conductivity reconstruction algorithm for MAT-MI, numerical simulations were conducted for a 2-layer eccentric cylindrical model with the conductivity similar to the conductivities of biological tissues. For simplicity, the acoustic wave was assumed to propagate in acoustic homogeneous media without considering the acoustic attenuation or reflection or scattering, and the sound speed is set to be 1500 m/s in all media. The cross-sectional conductivity distribution of the cylindrical model is shown in Fig. 2. The outer and inner media were set to have radii of 30 mm and 18 mm, centered respectively at (0, 0) mm and (−6, 0) mm where the corresponding conductivities are 1 S/m and 0 S/m. The scanning radius of the transducer was set to be 40 mm centered at (0, 0) mm. The entire system was placed in distilled water with the electrical conductivity of 0 S/m.
The acoustic waveforms collected by the transducers placed around the cylindrical model at the angles ranging from 1° to 360° were simulated. As indicated in Fig. 2, the distribution of acoustic pressure detected at (−40, 0) mm is illustrated in Fig. 3(a), which is normalized by the maximum pressure. Two positive pressure peaks (at 
and 
) and two negative pressure peaks (at 
and 
) are clearly displayed, representing the four tissue boundaries (A, B, C, and D in Fig. 2) with opposite conductivity variation directions. Due to the influences of the reception angle (
of the transducer and the radiation angle (
of the acoustic dipole sources, continuous pressure attenuations are generated, forming gradually pressure varying from one peak to zero and then to the opposite one. The positive and negative peaks A and D are generated by the dipole sources respectively at the boundaries A and D of the outer medium, while the peaks B and C are produced by the boundaries of the inner medium. The positive polarity of peaks (A and C) represents the conductivity change from 1 S/s to 0 S/m along the transmission path, while the negative ones (B and D) denote the conductivity change from 0 S/m to 1 S/m. The time interval between the peaks A and D is 
, which equals the diameter 60 mm of the outer medium. The vertical pressure jumps at the peaks A, B, C, and D exactly denote the BNPs produced at the corresponding boundaries, comprising the intensities and polarities of the vibration sources. In addition, although the strengths of the acoustic sources at the boundaries A and D are the same, the amplitude of peak A is higher than that of peak D because of the shorter transmission distance.
After convolution between the detected acoustic pressure and the system transfer function, the collected acoustic waveform was simulated and plotted in Fig. 3(b), which is normalized by its maximum amplitude. Corresponding to the four pressure peaks in Fig. 3(a), four wave clusters are clearly displayed. Meanwhile, the wave clusters generated by the continuous pressure attenuations are counteracted to zero. The amplitudes and vibration polarities of wave clusters are generated by the corresponding vertical pressure jumps (BNPs), reflecting the values and polarities of conductivity changes at the corresponding boundaries. The clusters A and D are in opposite vibration polarities with a time interval of 
, which denote the two boundaries between the outer medium and surrounding water. The clusters A and C produced by the corresponding positive pressure peaks are in the same vibration polarity, which is opposite to the vibration polarity of clusters B and D generated by negative pressure peaks.
After applying Hilbert transform to the collected waveform in Fig. 3(b), the envelope distribution of which was calculated and plotted in Fig. 4(a). Four envelope peaks A, B, C, and D produced by the wave clusters corresponding to the four conductivity boundaries are clearly displayed. With the analyses of amplitude and phase for wave clusters, four BNPs with corresponding amplitudes and polarities were retrieved, and the distribution of BNP along the measurement direction was restored as plotted in Fig. 4(b).
In order to realize the performance comparison with theoretical results, differential calculation for the simulated acoustic pressure 
in Fig. 3(a) was also performed and the normalized pressure derivative 
was calculated as illustrated in Fig. 5(a). Four dual-polarity spikes are located at the conductivity boundaries with “+/−” and “−/+” presenting the positive and negative polarities, respectively. Meanwhile, the absolute amplitudes of vertical pressure jumps were used to retrieve the corresponding BNPs. Through the analyses of amplitude and polarity for Fig. 5(a), the theoretical distribution of 
was restored as plotted in Fig. 5(b). It is obvious that the positions, amplitudes and polarities of vertical pressure jumps (BNPs) for spikes A, B, C, and D in Fig. 5(b) are basically the same as those in Fig. 4(b), demonstrating that the proposed algorithm using Hilbert transform can be used to restore the distribution of BNP from the collected waveform successfully as a simplified method.
With the 360 simulated waveforms collected around, the traditional MAT-MI image was reconstructed as shown in Fig. 6(a). Two eccentric rings with the radii of 30 nm and 18 mm are displayed, which are consistent with the cross-sectional conductivity distribution in Fig. 2 in terms of shape and size. Two fringe boundaries in opposite polarities (inner: −/+/−, outer: +/−/+) reflect the reverse vibration polarities of the wave clusters produced by the boundaries with opposite conductivity changes (0 S/m 
1 S/m and 1 S/m 
0 S/m) from inside to outside. The spatial resolution of the reconstructed image is reduced obviously by the wide width (
mm) of the boundary fringes, which are determined by the system transfer function. However, besides the boundaries, the conductivity distribution inside the media cannot be visualized in Fig. 6(a).
By applying the proposed reconstruction algorithm to the collected waveforms, 360 distributions of BNP were restored. Based on Eq. (13), the 1D relative distribution of 
at each measurement angle was calculated and the 2D image of normalized conductivity was reconstructed by using Eq. (14) as plotted in Fig. 6(b). Besides the shape and size, the conductivity distributions inside the inner and outer media show good agreement with those in the model. Whereas, obvious conductivity difference inside each medium can still be visualized on different gray scales, which is introduced by the superposition effects[19,28] during the tomographic calculation in image reconstruction.
4. DiscussionAs is well known, the acoustic source in MAT-MI is generated by the derivative of the exerted Lorentz force, which is determined by the conductivity gradient along the normal direction of the induce eddy current.[31] Higher acoustic pressure can be generated at the boundary with sharper conductivity variation. Therefore, in most of previous studies, cylindrical phantom models with sharp conductivity boundary are often chosen to calculate source distribution for the scanned layer. Only the effects of boundary sources are taken into account without considering the inner sources. With the simulations as mentioned above, it can be concluded that the proposed BNP based reconstruction algorithm can be used to reconstruct the conductivity distribution of the scanned layer in terms of shape and size. And comparing with traditional MAT-MI, the conductivity distribution inside the media can also be reconstructed with improved spatial resolution and enhanced image definition. Whereas, for the practical measurement of early-stage or invasive cancer, conductivity gradual-varying distribution is proved to be existent at the boundary between the normal and pathological tissues.
In previous studies of our laboratory, acoustic source analysis for MAT-MI is conducted for biological tissues with conductivity gradual-varying boundaries.[31] It is proved that the inner source is generated by the product of the conductivity and the curl of the induced electrical intensity, while the boundary source is produced by the cross product of the gradient of conductivity distribution and the induced electrical intensity. For biological tissues with low conductivity, boundary source is determined by the conductivity variation rate, while the inner source is determined by the induced electrical field. Only for sharp conductivity boundaries, will the strength of boundary source be much higher than that of the inner source and the influence of inner source can be neglected. Therefore, for biological tissues with conductivity gradual-varying boundaries, BNPs are too low to be retrieved from the low-level wave clusters, when only applying the proposed simplified BNP based conductivity reconstruction algorithm. Further characteristic analyses for the conductivity gradual-varying model exhibit great significance for the practical application of MAT-MI, and new reconstruction method with more efficiency and precision should be developed in future studies.
In addition, in order to analyze the spatial resolution of the proposed conductivity reconstruction algorithm, several simulations are also conducted for the cylindrical model by left moving the position of the inner medium. It is proved that with the distance between the boundaries A and B decreasing from 6 mm to 1 mm, the clusters B and C move left, and obvious cluster superimposition between the clusters A and B can be observed to generate a wider wave cluster with relatively enhanced amplitude when the distance is less than 3.0 mm (wave length). For the cross-sectional conductivity distribution of the model with a distance of 2 mm between the boundaries A and B as shown in Fig. 7(a), the 2D conductivity image of the scanned layer obtained with the proposed reconstruction algorithm is illustrated in Fig. 7(b). The reconstructed image shows high consistency with the cross-sectional distribution of the model. Whereas, as the distance between the boundaries A and B is less than 3 mm, obvious errors can be visualized with aliasing boundaries (A and B) and artifacts in surrounding water. Therefore, although BNPs can be retrieved from the envelope peaks and vibration phases of the collected waveforms by using Hilbert transform, is the spatial resolution of the proposed simplified BNP-based conductivity reconstruction algorithm demonstrated to be one wavelength due to the superimposition influence of wave clusters.
5. ConclusionsIn this paper, based on the analyses of MAT-MI, the inverse problem solution for conductivity reconstruction is derived in explicit formula with clear physical meanings of acoustic pressure and the corresponding pressure derivative. A simplified BNP based conductivity reconstruction algorithm is proposed for the collected waveforms by using Hilbert transform. Numerical simulations of acoustic pressure and conductivity reconstruction are conducted for a 2-layer eccentric cylindrical model, and also compared with the results of traditional MAT-MI. The good agreement between the reconstructed images and the cross-sectional distributions of the model demonstrates the feasibility of accurate conductivity reconstruction with improved spatial resolution (one wavelength) and image contrast. This study provides a simplified method of accurate conductivity reconstruction for MAT-MI and suggests potential applications in imaging pathological states of tissues with conductivity variations.
